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ABSTRACT 

Acceleration and transport of high-energy particles and fluid dynamics of atmospheric plasma are 
interrelated aspects of solar flares, but for convenience and simplicity they were artificially separated in 
the past. We present here self-consistently combined Fokker-Planck modeling of particles and hydro- 
dynamic simulation of flare plasma. Energetic electrons are modeled with the Stanford unified code of 
acceleration, transport, and radiation, while plasma is modeled with the Naval Research Laboratory 
flux tube code. We calculated the coUisional heating rate directly from the particle transport code, 
which is more accurate tha n those in prev i ous st udies based on approximate analytical solutions. We 
repeated the simulation of IMariska et al.l ()I989l ) with an injection of power-law, downward beamed 
electrons using the new heating rate. For this case, a ~10% difference was found from their old result. 
We also used a realistic spectrum of injected electrons provided by the stochastic acceleration model, 
which has a smooth transition from a quasi-thermal background at low energies to a nonthermal 
tail at high energies. The inclusion of low-energy electrons results in relatively more heating in the 
corona (vs. chromosphere) and thus a larger downward heat conduction flux. The interplay of electron 
heating, conduction, and radiative loss leads to stronger chromospheric evaporation than obtained in 
previous studies, which had a deficit in low-energy electrons due to an arbitrarily assumed low-energy 
cutoff. The energy and spatial distributions of energetic electrons and bremsstrahlung photons bear 
signatures of the changing density distribution caused by chromospheric evaporation. In particular, 
the density jump at the evaporation front gives rise to enhanced emission, which, in principle, can 
be imaged by X-ray telescopes. This model can be applied to investigate a variety of high-energy 
processes in solar, space, and astrophysical plasmas. 

Subject headings: acceleration of particles — hydrodynamics — methods: numerical — Sun: 
chromosphere — Sun: flares — Sun: X-rays, gamma rays 



1. INTRODUCTION 

A solar flare, as one of the most prominent manifes- 
tations of solar activity, has many faces among which 
are acceleration and transport of high-energy particles 
and the dynamical response of atmospheric plasma. It 
is generally believed that magnetic reconncction in the 
corona is the primary energy release mechanism that 
leads to plasma heating and particle acceleration. The 
heated plasma and accelerated particles (primarily elec- 
trons) produce bremsstrahlung X-rays at the apex of 
the flare loop observed as a loop-top (LT) source (e.g., 



Masuda et ai.lll994HPetrosian et al.ll2002l:lKrucker fc LinI 



20081) . Some of the released energy is transported down 
the closed magnetic loop by nonthermal particles (elec- 
trons and ions) and thermal conduction, which con- 
tribute to energy gain in various layers of the atmosphere. 
Electrons give up most of their energy to ambient par- 
ticles via Coulomb collisions, and produce hard X-rays 
(HXRs) primarily at the footpoints (FPs) of the loop in 
the dense transition region (TR) and chromospher e (see 
lHovng|[l98ll : lSakaol[T99i : iSaint-Hilaire et al.ll200a[ ). Ac- 
celerated protons and heavier ions, on the other hand, 
cause nuclear reactions while colliding with background 
particles and produce 7-rays. Some accelerated electrons 
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and ions escape along open magnetic field lines into the 
interplanetary space and are observed as sola r energetic 
particles (SEPs) by in situ instrum ents (e.g.. lLinllI985l : 
iLiu et al.ll2004atlKrucker et al.ll2007D . The rate of energy 
gain in the chromosphere, if exceeding the combined lo- 
cal radiative and conductive cooling rate, can rapidly 
heat the plasma up to a temperature of ^10'' K. The 
resulting overpressure drives an upward mass flow at a 
speed up to hundreds of km s~^, which fills the flare loop 
with a hot plasma, giving rise to the gradual increase of 
soft X-ray (SXR) emissio n. This process , termed chro- 
mospheric evaporation bv lNeuperd (|1968f l. can influence 
particle transport by changing the ambient density in 
the loop on timescales of tens of seconds, and affect heat 
conduction by changing the loop temperature distribu- 
tion at the same time. CoUisional and conductive heat- 
ing will be consequently modified, and in turn, so will 
the dynamic atmospheric response. On longer timescales 
of minutes, as magnetic reconncction proceeds and new 
loops are formed or excited, the above processes repeat 
sequentially in newer loops. 

1.1. Motivation for This Study 

The aforementioned processes are coupled in a circu- 
lar chain, but due to great complexity of the subject, 
previous researchers tended to focus on one process at a 
time while assuming some simple forms for others. Past 
investigations fall into two general categories: (1) accel- 
eration and transport of particles and (2) fluid dynamics 
of atmospheric plasma. 
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(1) Various mechanisms have been proposed for 
particle acceleration. Among the agents o f ac- 
celeration are DC electric fields (iHolmanI 119851: 
lLitvinenkolll996t IZharkova fc Gordovskvvll2004f) . shocks 
(iTsuncta & Naito"1998')^ and turbulence or plasma waves 
(IRamatv.. 1979- Hamilton fc Petrosianlll992t iMiller et all 
119961 iPetrosian fc Liull2004f) . Particle transport IS com- 
paratively better understood and previous studies usu- 
ally assume d a hydrostatic atmosphere. Early analyti- 
cal st udies (Brown 1973': P etrosianl[l97a iLin fc HudsonI 
119761: fEmslic 197 8: Ch audrashekar fc Emslidll987t) took 
various approximations (e.g., neglecting pitch-angle dif- 
fusion) to allow the problem to be tractable. The nu- 
merical study of iLeach fc PetrosianI ()1981[ ) improved on 
this by solving the Fokkcr-Planck transport equation 
with inclusion of pitch-angle changes due to Coulomb 
collision and magnetic mirroring. This was later ex- 
tended to the relativistic regime by including energy 
losses a nd pitch-angle changes due t o synchrotron ra- 
diation (|McTiernan fc Petrosianlll99Cl[ ). Similar Fokker- 
Planck studies of pa r ticle t ra nsport were perform ed by 
MacKinnon fc Craid (1199 iD iMcClementJ (|1992[ ). and 



Svniavskii fc Zharkoval ()1994[ ) 



(2) Fluid dynamics of the magnetized atmosphere 
in response to flare heating can be best studied 
with a multi-dimensional magnetohydrodynamic (MHD) 
model, but for simplicity most efforts were invested in 
one-dimensional hydrodynamic (HD) simulations. This 
is justified for the solar corona where the magnetic pres- 
sure dominates the gas pressure (low (3 plasma) and 
resistivity is low. As a result, plasma is only allowed 
to flow along the magnetic field lines due to the line- 
tying c ondition. Previous HD mode ls ('Mac Neice et al.l 
1981 INaeai fc Emslid [l98i: lEl^ lic & N a"gail [l985l: 



Fisher ct al."1985a(lb'!d:l Mariska et arill989[ iGan fc Fand 
1990: Emslic ct a l. .1993) usually assumed a power-law 
spectrum of accelerated electrons injected at the apex 
of the loop, and calculated coUisional heating along 
the loop by these electrons from approximate analyt- 
ical solutions of p articl e tran sport mentioned a bove. 
lAbbett fc Hawlevi (fT999h and lAUred et all pOOSh im- 
proved on previous studies by including detailed calcula- 
tion of radiative transfer in the atmosphere. 

There are theoretical and observational motivations to 
investigate the particle and fluid aspects of a solar flare 
together in a self-consistent manner. From a theoretical 
point of view, such an investigation is demanded in or- 
der to retrieve missing physics when the two aspects were 
studied separately. It has also become technically more 
feasible, thanks to advances in both aspects over the 
last three de cades and particularly in recent years. Sev- 
eral studies (jMiller fc Mariskall2005l : I Winter et alll2007f) 
along this direction are already under way, but none of 
them has been completed. From an observational point 
of view, new observations, particularly X-ray images and 
spectra obtained by the current Ramaty High Energy So- 
lar Spectroscopic Imager (RHESSI) and previous Yohkoh 
missions, have posed new challenges to t he existing thco - 
ries. For example, in recent stud ies of thelNeupertI (|1968f ) 
effect. IVeronig et al.l (I2005D and lLiu et all (|2006bD found 
that, unexpectedly, the electron energy deposition power, 
which is more physically related to the thermal energy 
change rate, did not yield a better correlation with the 
time derivative of the SXR flux than the conventional 



HXR flux. In an event of chromosphe ric evapora t ion im - 
aged by RHESSI for the first time, iLiu et al.l (|2006bD 
found X-ray sources moving from the FPs to the LT 
at very h i gh spe eds (~10^ km s~^). More interestingly, 
ISui et ail ()2006l ) found double nonthermal sources mov- 
ing first downward from the LT toward the FPs and then 
upward along the loop. To fully understand these obser- 
vations requires a joint study of acceleration and trans- 
port of particles and fluid dynamics of the atmospheric 
response. 

1.2. Approach of This Study 

With the goal to investigate the coupled pro- 
cesses of acceleration, transport, and hydrodynam- 
ics in solar flares, we present here combined Fokkcr- 
Planck modeling of particles and HD simulation of 
plasma. (1) The Fokker-Planck model utilizes the 
Stanford unified code of particle acce leration, trans- 
port, and bremsstrahlung radiation (jPetrosian et al.l 
|2001| ) . The transpo rt and radiation calcula t ion is 
base d on the work of iLeach fc PetrosianI (jl98ll I1983D 
and iMcTiernan fc PetrosianI "( 199(3[ ). The acceleration 
module of the code adopts the stochastic acceleration 
model of IPetrosian fc Liul (|2004l hereafter PL04), which 
has inherited k nowled ge ac c umulated over a decade 
(iHamilton fc Petrosian. 1992t iDung fc PetrosianI 1 19941 : 
iPark fc PetrosianI 119951 . il996l : iPark et al.lll997t l. When 
compared with observations, th is model has many at- 
tract ive features and advantages (|Liu et al.ll2004al . r2006al . 
l2008l) over other mechanisms. (2) The HD simula- 
tion uses the Naval Research Laboratory (NRL ) So- 
lar Flux Tube Model (jMariska. Emshe. fc Lilll989l . here- 
after MEL89 ) , whi ch, as a modified version of the 
iMariska et al.l ()1982l ) model, provides excellent treat- 
ment of fluid dynamics and has been widely used in 
studying atmospheric re sponse to flare heating (e.g., 
iWarren fc Doschekll2005l ). 

One of the major advances marked by this study is 
the more accurate and self-consistent evaluation of the 
collisional heating rate by nonthermal electrons. This 
heating rate is critical to HD simulation of flares, but 
was not properly calculated previously in two major as- 
pects: (1) The calculation of energy loss of energetic 
electrons and thus the heating rat e was based on ap- 
proxi mate analytical solutions (e.g.. lBrownlll973tlEmslid 
|1978() . which incorporated only pitch angle growth due 
to Coulomb collisions, but in reality the pitch angle 
change is a diffusion process. This will be remedied 
in this study with the inclusion of a full Fokker-Planck 
treatment of electron transport. (2) Another previous 
drawback was the use of an unrealistic spectrum of in- 
jected electrons, which usually was a power-law with a 
low-energy cutoff. iFisher et al.l ()1985cD . for example, as- 
sumed a sharp low-energy cutoff at Ec = 20keV (i.e., no 
electrons below i?c), while MEL89 introduced a "soft" 
cutoff below which the spectrum is still a power-law 
with a positive slope. It should be noted that rather 
than an intrinsic property of the primary accelerated 
electron population, a low-energy cutoff or turnover in 
the electron spectrum inferred from X-ray observations 
can re sult from secondary effects su ch as return cur- 
rents (IZharkova fc Gordovskvvl l2006f) and photospheric 



albedo (iLanger fc PetrosianI Il977t iBai fc Ramatvl Il97"8l : 
ISui et al.ll20'07ir The collisional heating rate is sensitive 
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to the injected electron spectrum and thus the use of an 
incorrect spectrum would make the HD simulation devi- 
ate from reality significantly. PL04 has provided a more 
realistic electron spectrum that has a continuous span 
from a quasi-thermal distribution at low energies to a 
nonthermal tail at high energies, avoiding unnecessary 
low-energy cutoff. It also gives good fits to both LT and 
FP X-ray spectra obtained by RHESSI. Such an electron 
spectrum is used in this work. As we will see later, the 
low-energy electrons, which otherwise would have been 
missing if a cutoff were to be present, play an impor- 
tant role in heating and in influencing the subsequent 
HD evolution. 

We present the numerical model in § [21 techniques to 
combine different modules of the model in § |3l and simu- 
lation runs in § S) We compare the HD characteristics of 
different simulations in § [5] and examine the HD effects 
on particle transport and X-ray emission in § [6l Con- 
clusions and discussions are given in [71 This model is 
based on the PhD t hesis ofjLiul ()2006l ). which has also 
appeared as a book ()Liull2008[ ). The model has been re- 
fined ever since and the numerical results presented here 
are new. Our new calculations and analyses include: (1) 
recomputing the MEL89 simulation using our transport 
code (see ij4.2[) . (2) detailed analysis of the interplay of 
heating and cooling (i i5.2p . and (3) examining the tem- 
perature distribution of plasma velocity which can be 
directly compared with Doppler observations ( iJ5.3p . In 
order to achieve more accuracy, we have extended the 
number of pitch angle bins from 24 to 100. 

2. SIMULATION MODEL 

Here we consider the dynamical evolution of a single 
flare loop perpendicular to the solar surface. As shown in 
Figure[Tl(top), an acceleration rej^ion of length L = 5 Mm 
is located at the top of the model loop, sandwiched be- 
tween two symmetric quarter-circles called the transport 
region of length Smax = 14 Mm. The loop has a uniform 
circular cross-section a(.s) = ao = tt?'^ with a constant ra- 
dius of r = 0.3 Mm at any distance s from the edge of the 
acceleration region. In its initial state (Fig. [H bottom) 
the loop spans from the hot [T > 10^ K), tenuous corona 
to the cold (T = 10^ K), dense chromosphere, with the 
TR (defined here as the lowest point where T > 10^ K) 
located at str — 10 Mm. 

The simulation model includes four parts: (1) The 
stochastic acceleration code generates a spatially aver- 
aged distribution (in energy) of high-energy electrons in 
the acceleration region. This distribution is fed to the 
transport region where (2) the transport code computes 
the electron distribution (in energy and pitch angle) and 
collisional heating rate as a function of distance (depth), 
and (3) the hydrodynamic code simulates the atmospheric 
response to this heating. Finally, (4) the radiation code 
calculates the corresponding bremsstrahlung emission in 
the loop. Parts (1), (2). and ( 4) are inherited from the 
Stanford unified particle code (jPetrosian et al.ll2001[ ) in 
which the acceleration module has been revised accord- 
ing to PL04, while part (3) is adopted from the NRL 
flux tube code (MEL89). Details of the four parts are 
described in the following subsections. 

Note that sequential energizing (or excitation) propa- 
gating from one flare loop to another h as been observed 
as the apparent motions of HXR LT ([Gallagher et al.l 
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Fig. 1. — Top: Geometry of the model flare loop. 8 is the pitch 
angle of the electron with a velocity Ve in the guiding magnetic 
field Bo- Bottom: Initial distribution of logarithmic temperature 
T, gas pressure P {left scale), and electron number density rie (right 
scale) vs. distance in the transport region of the loop, f'ressure is 
scaled upward by a, factor of 100. 

200alSui fc Holmanll200l a nd FP (iGrigis fc Benzl[2005l : 
Yang et all [20091: iLiu et^ [20091) sources. In our sim- 
ulations, we set the duration of the impulsive p hase to 
be 60 s. This time, according to [Schriiver et all (|2(306il . 
translates into an apparent FP speed of = lOkms"^ 
(where r = 0.3 Mm), which is comparable to the ob- 
served values ([Liu et al.[[2004br ). Our single loop scenario 
is thus legitimate within the relevant timescale and can 
be viewed as an elementary process of sequential exci- 
tation of multiple loops. Evolution on longer timescales 
(say, >100 s) involves multiple loops and can be stud- 
ied b y superimposi ng sequential single-loop simulations 
(e.g.. [Warre'^ [2001 . 

2.1. Stochastic Acceleration 

The stochastic acceleration model of PL04 addresses 
electron and proton acceleration by plasma waves prop- 
agating parallel to the background magnetic field Bq. 
According to this model, large-scale turbulence or long- 
wavelength plasma waves are generated in the corona as 
a result of magnetic reconnection. The turbulence, cas- 
cading to smaller scales, heats plasma and accelerates 
particles in a region near the top of the flaring loop. The 
heated plasma and accelerated electrons produce the ob- 
served thermal and nonthermal X-rays, respectively, in 
the acc eleration region or the LT source (IXu et al.ll2008l : 
[Liu et a l. 2008). Here we briefly repeat the mathemat- 
ical description of the model. The Fokker-Planck equa- 
tion that governs electron acceleration, in general, can 
be written as 



dt 



d_ 

dE 



D{E) 



dE 



d_ 

'dE 



{[A[E)-El]U} 
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Tcsc{E) 



Q{E), 



(1) 
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where /ac = fadt^E), in units of electrons cm^'^keV"^, 
is the angle-integrated and spatially averaged electron 
distribution function in energy space (subscript "ac" de- 
notes the acceleration region), E is the electron kinetic 
energy, D{E) is the energy diffusion coefficient, A{E) the 
direct acceleration rate, TcsciE) the particle escape time, 
Q{E) is the rate of background electrons being supplied 
to the acceleration region that serves as a source term, 
and finally El = -Ecoui + £'synch is the absolute value of 
the net energy loss rate that is a sum of the Coulomb 
and synchrotron loss rates. 

In order to solve equation ([T]) one must first eval- 
uate all the terms. The energy loss rates E'coui and 
-E^synch are well known (see eq. [18] of PL04) and the 
source term Q{E) is to be prescribed with a specific 
model (assumed to be a thermal or Maxwellian distri- 
bution here). The central task is then to obtain D{E), 
A{E), and Tcsc{E). They are determined through the 
dispersion relation of the plasma waves (PL04 eq. [28]) 
and the wave-particle resonance condition (PL04 eq. [4] ) . 
Following PL04, we assume a fully ionized H and ^He 
plasma with a relative abundance of electron/proton/a- 
particle=l/0.84/0.08, and a broken power-law spectrum 
of the turbulence given by equation (29) in PL04 with the 
relevant spectral indices qi = 2, q = —1.7 (Kolmogorov 
value), and q^ = —4. The characteristic acceleration rate 
given by equation (30) in PL04 represents the rate 
of wave-particle interaction and depends on the level of 
turbulence. 

Once all the coefficients have been evaluated, equation 
p]) is solved numeric ally using the flux conserv ative finite 
di fference scheme of IChang fc Coopeij ()1970() described 
in iPark fc PetrosianI ()1996f ). Here we assume a homoge- 
neous acceleration region and obtain a steady state solu- 
tion of fac{E) (i.e., d/dt = 0). The angle-integrated flux 
in the acceleration region is Fac{E) — Vefac{E), where 
f e is the electron velocity. We then calculate the flux of 
electrons that escape the acceleration region and enter 
the transport region of the flare loop. 



Fcsc{E) 



UcjE) 
TcsciE) 



(2) 



2.2. Particle Transport 



The flux Fcf^r.(E) is then input to the particle transport 
code (iLeach fc PetrosianI Il981l : iMcTiernan fc PetrosianI 
|199C1() which calculates the electron distribution in en- 
ergy and pitch-angle space and its variation with dis- 
tance while the electrons spiral down magnetic field lines 
into deeper layers of the atmosphere. The code numer- 
ically solves the fully relativist ic, steady-state, Fokker- 
Planck equation (i.e., eq. [1] in IMcTiernan fc PetrosianI 
|199CI() . which is similar to equation ^ here and includes 
energy loss (no energy diffusion) and pitch angle diffu- 
sion due to Coulomb collision, and pitch angle changes 
due to magnetic mirroring and synchrotron radiation. 
Following ^ IcTic man ( 198JD, we neglect return cur rents 
(|Svniavskii fc Zharkovalll9M iZharkova et"aI1[l995l) . 

The variable^ to be solved in the transport equation 
is the electron flux spectrum F(E^ fi, s) as a function of 

In practice, the numerical code equivalently solves 
for F{E,^,s)/l3^ = c$a(s)/ao, where <& is defined in 
IMcTiernan fc PetrosTaiil II1990I ) and 13 = Ve/c. 



energy E, cosine fj, — cos 6 of pitch angle 0, and dis- 
tance s from the injection site at the boundary of the 
acceleration region. F{E, fi, s)dfj, has units of electrons 
s~^ cm~^ keV~^ and is evaluated as 

F{E,^i,s)^vJiE,^,,s)^, (3) 

ao 

where f{E, ^, s)dfi is the number density distribution 
function in units of electrons cm~'^ keV~^ (cf., the angle- 
integrated number density f^dE) in the above accelera- 
tion code) , and we integrate the differential electron flux 
Vef(E, ^, s) over the cross-sectional area a(s) of the loop 
and then normalize it by a constant equivalent area oq. 
As noted earlier, here we assume a constant a(s) = oq for 
simplicity, which means a uniform magnetic field along 
the loop and thus no magnetic mirroring. 

In addition to the injected electron flux FcsciE) from 
the acceleration code, the transport code requires the 
knowledge of the ambient density and abundance dis- 
tribution along the loop. (1) Here we assume that 
i^csc is isotropic in pitch angle, representing the con- 
sequence of frequent scatterings of electrons by tur- 
bulence in the acceleration region. This assumption 
is consistent with the nearly isotropic, rather than 
beamed, distributions inferred from center-to-limb vari- 
ations of HXR and 7-ray fluxes and spectral indices in 
observations obtained by the Solar Maximum Mission 
(jMcTicr nan & Pctrosian 1991.) , and more recently from 
atmospheric al bedo due to Compton back-scattering in 
RHE SSI flares (jKontar fc Brownll2006l : iKasparova et all 
|2007| ). The injected flux at the top (s = 0) of each leg of 

the loop is then F{E,h,s)\s=q = Fesc{E)/ J^-^dfi which 
is equivalent to a uniform distribution in the Air solid an- 
gle integrated over the 27r range of the azimuthal angle (/) 
assuming axis-symmetry. With the symmetric assump- 
tion, the steady state calculation is performed in only one 
leg of the loop. We impose a symmetric (or reflective) 
boundary condition at s = 0, where a particle leaving 
the computational domain is reflected back with identi- 
cal energy but opposite pitch angle cosine, mimicking a 
particle coming from the other leg of the loop. (2) As to 
the background atmosphere, we assume a fully ionized 
hydrogen plasma whose distribution is taken from the 
result of the HD code described next. 

2.3. Hydrodynamics 

Hydrodynamics in the transport region is calculated 
with the NRL solar flux tube code (MEL89) based on 
iMariska et all (|1982[ ). The code assumes a two-fluid 
plasma composed of electrons and ions that can only 
move along the magnetic field in a flux tube. The user- 
speciflcd geometry of the tube (a uniform quarter-circle 
in our case) is characterized by the tube cross-sectional 
area a(s) and the component of the gravitational ac- 
celeration along the tube. The code solves the time- 
dependent, one-dimensional equations of mass, momen- 
tum, and energy conservation (see eqs.[l]-[3] in MEL89). 
The independent variables are the mass density p(s), 
fluid velocity v{s), and temperature T(s) which we as- 
sume to be the same for electrons and ions. Because of 
small masses of electrons, we neglect the momentum loss 
of energetic electrons to the background plasma. 

The volumetric heating rate in the energy equation is 
S{s)^Sds) + S^, (4) 
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where Se{s) represents heating by energetic electrons, 
which is provided by the transport code (see |J3|), and 
So = 8.31 X lO-^ergss"^ cm'^ (MEL89) represents uni- 
form background heating, presumably caused by coronal 
heating in the quiet sun active region. The conductive 
fl-ux -Fcond and heating rate Scond arc 
dT 



cond 



s, 



cond 



(5) 



ds ' 
The radiative en- 

(6) 



where k is the thermal conductivity, 
ergy loss rate is 

where rie and Up are the electron and proton number den- 
sity, respectively, which are equal by our assumption of 
fully ionized hydrogen plasma, and <I>(T) is the optically- 
thin radiative loss function (MEL89) which has its max- 
imum at T ~ 1 - 3 X lO^K. 

We select an adaptive mesh of 450 grids that move with 
time to optimize spatial resolution in the dense chromo- 
sphere and near sharp jumps at the TR and evaporation 
front. This mesh is also shared by the transport and 
bremsstrahlung radiation codes in our new model. A re- 
flective (or symmetric) boundary condition is imposed 
at both the upper (loop apex) and lower (deep in the 
chromosphere) boundaries of the transport region (see 
Fig.[T]), such that the system remains closed. 

2.4. Bremsstrahlung Radiation 

Having obtained the electron flux from the transport 
code and the background density from the HD code, 
we calculate the thin-target bremsstrahlung radiation 
intensity or photon emission rate, /(e,s), as a func- 
tion of photon energy e and distance s. /(e, s) (pho- 
tons s~^ em"-^ keV^^) is defined as 



I{e,s) 



dE 



np{s) 



d(T{e, E) 
Je 



(7) 



where da{e,E)/de is the angle-averag ed differential 
bremsstrahlung cross-section given by iKoch fc Mot j 

(|1959f) . and Fint = J^-^^ F{E, s, fi)dfi is the angle- 
integrated electron flux. Substituting Fint with the accel- 
eration region flux i^ac gives the LT emission /LT(e), while 
identifying Fmt = LFthick yields the spa.tially integrated 
thick -target emission /thick(e) (jBrownl Il971l iPetrosiaiil 
Il973l ). Here -Fthick = "f^e/thick is the equivalent thick- 
target el ectron flux, given by the co rresponding number 
density ijPetrosian fc Donaghvl[T999l : PL04) 
1 f°° 

fthick{E) = -^ F,,,{E')dE' . (8) 
LEl je 

The equivalent FP emission is the spatially averaged pho- 
tons emitted below the TR (located at St, ), /Fp(e) = 
Xf™ /(e, s)ds/(sniax - Str), whcrc Sniax IS the distance 
at the lower boundary of the loop. If the coronal is neg- 
ligibly tenuous and the column depth at Smax is large 
enough to stop all HXR producing electrons of interest, 
sti)lFp{£) approaches /thick(e)- 
Comparison between the HXRs modeled here and 
those observed by RHESSI satellites can serve as a 
unique diagnostic tool and will be pursued in a future 
publication. Here we show an example (Fig. [2]) of how 
well observed LT and FP fluxes can be fitted with the 
above equations using a spectrum of accelerated electrons 
given by the PL04 model. 
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Fig. 2. — Photon fluxes multiplied by energy squared (e^) during 
the main HXR peak of the 2002 August 3 XI. flare and spectral 
flts using the stoehastic acceleration model of PL04. ( a) Summe d 
fluxes of all LT sources and all FP sources (see Fig. 2.6 of lLiul2006l '). 
(6) Sum (squares) of the LT and FP fluxes shown in (a) and sum 
{solid line) of the corresponding model fits. Overplotted are the 
spatially integrated spectrum [plus signs) and the corresponding 
preflare background [dotted line). The legend lists parameters used 
in the model. (We thank Siming Liu and Yan-Wei Jiang for help 
in producing this figure.) 

3. COMBINING PARTICLE AND HYDRODYNAMIC CODES 

The main task in this study is to combine the Stanford 
particle code and the NRL HD code. Here we assume 
that particle acceleration acts as an independent driver of 
the simulation and is not affected by the temperature or 
density evolution in the transport region of the loop. The 
task is thus reduced to make the transport module of the 
particle code and the HD code communicate interactively 
in real time. 

3.1. Electron Heating Rate 

The rate of coUisional heating to the background 
plasma, Se (ergs s~^ em~^), equals the rate of energy 
loss from the energetic electrons. This can be calcu- 
lated from the electron distribution given by the trans- 
port code in two equivalent ways, the second of which is 
used in this study. 

(1) Se can be evaluated from the energy loss rate i?coui 
due to Coulomb collisions as 



Seis) = 



dE 



f{E, iJ.,s)Ecou\dfi , (9) 



where [i^min, F'max] is the range of the energy bins used 
in the simulation, and the electron distribution function 
f{E,fi,s) can be obtained from the corresponding elec- 
tron flux F{E,fi, s) via equation ([3]). 

(2) Alternatively, one can calculate the net downward 
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energy flux carried by the electrons, 

Fcrg{s) = -^ [ dE I ^iEF{E,s,^i)d^i, (10) 

J-i 

and differentiate it to obtain the net energy gain to the 
background plasma in a unit volume, 

Seis)^~dF,,^is)/ds, (11) 

where iiEF{E, s, /i) is the energy flux along the loop and 
the factor ao/a{s) accounts for the variation of the loop 
cross-sectional area. This approach is practically equiv- 
alent to equation ([9]), because in the HXR energy range 
the combination of synchrotron and bremsstrahlung ra- 
diation only constitutes a negligible fraction (< 10""') of 
the energy loss of a fast electron due to Coulomb colli- 
sions. 

3.2. Code Communication 

It is desirable that the particle and HD codes com- 
municate at each time step during the time advance. 
The current transport code, however, can only provide a 
steady-state solution and does not have a time-dependent 
capability. This can be remedied by carefully select- 
ing the communication time interval A^, because par- 
ticle transport occurs on a much shorter timescale than 
hydrodynamics. This interval should be as short as pos- 
sible provided that a steady-state transport solution can 
be reach ed. By co nsidering the electron "lifetime" (see 
eq. [9] of lPetrosianI 1973), which is determined by the en- 
ergy loss time in a gi ven l oop g eometry and atmospheric 
density distribution. iLiul (|2006l . his §7.2.4) found the op- 
timal interval to be At = 2 s. 

The remaining question is what heating rate the HD 
code should use during its time advance between adja- 
cent communications with the particle code. Let us first 
change the independent variable of Se from distance s to 
column depth N{s) = ne{s')ds', 

Se[N{s)]^Se{s)/ne{s), (12) 

noting Seis)ds = SeiN)dN and dN = Ueds. Here Se{N) 
is in units of ergs s^^. We have assumed a loop of uni- 
form cross section and thus no magnetic mirroring, and 
here we further neglect synchrotron loss -Egynch (valid for 
<1 MeV electrons). Under these assumptions, the elec- 
tron flux F[E,^i,s(N)] = F[E,ii,N(s)] is a function of 
column depth N ^ independent of distance s, and so docs 
the heating rate Se{N) calculated from F{E, ^, N). 

The HD response timescale is characterized by sound 
travel time, which is 84 s for a sound speed of 166 kms~^ 
at T = 10*^ K in a 14 Mm long loop here. Since At = 
2s ^ 84s, we assume that Se{N) is constant in time 
between adjacent code communications. During this At, 
the spatial distribution of the heating rate Se{s,t) varies 
with time merely according to the redistribution of den- 
sity and thus the variation of column depth, 

Se{s,t) = SeiN)ne{s,t). (13) 

In practice, at a given time t and distance s, we first 
identify its column depth 7V(s,t), which is then used to 
evaluate the heating rate Se{N), and then we apply the 
local density to obtain Se{s,t) by equation (fT3|) . 

Communications between the two codes are summa- 
rized in Figure [3l First, the HD code passes the initial 



density distribution to the particle code, which then runs 
its first steady-state calculation and returns the heating 
rate Se{N) as a function of column depth to the HD 
code. Next the HD code repeatedly converts Se{N) to 
S'e(s, t) as a function of distance at each time step using 
the latest density profile. Once the HD code advances a 
time interval of At = 2 s, it passes the updated density 
distribution back to the particle code, which starts the 
next cycle of iteration. 





Start 




Density ne(s) 






Particle 1st run 


HD initial state 






Heating rate Se{N) 








Density ne(s) 






Particle run 


HD run At=2 s 





Fig. 3. — Task flow chart for communications between the par- 
ticle and HD codes. 



4. SIMULATION RUNS 

We have performed three simulation runs (see Table [1]) 
to test the relative effects of different processes. (1) In 
the first run, which we refer to as the Old Model (abbre- 
viated by "O"), we assumed an injection of electrons of 
a power law at the LT, and evaluated the heating rate 
and HD response as in the MEL89 model. This model 
does not calculate particle transport properly. (2) In the 
second simulation, we still injected power law electrons 
and evaluated electron transport and heating along the 
loop using our transport code. We call this the Hybrid 
Model (or "H"). (3) Finally, we employed our most re- 
alistic model, where we evaluated from our acceleration 
code (PL04) the spectra of electrons at the LT accelera- 
tion site and those escaping the LT region. We also cal- 
culated transport and heating using our transport code. 
We call it the New Model (or "N"). In all cases, we as- 
sumed an identical initial HD state as shown in Figure [TJ 
and calculated the HD response using the MEL89 code. 
We assumed the dynamic or modulation profile of the 
number of injected electrons (power law for Runs O and 
H and Q for Run N) to be a triangular shape with a rise 
and fall to be 30 s each. Beyond this first 60 s of the 
impulsive phase, while the electron heating rate Se was 
set to zero, we continued the computation into the decay 
phase until t = 90 s. 

4.1. Run O: Old Model 

We first computed the HD response using the same 
heating function and almost identical control parameters 
as the "Reference Calculation" of MEL89. This m odel 
is based on the analytical solution of lEmslii (|1978[ ) that 
includes only energy loss and pitch angle growth of in- 
jected electrons. By the MEL89 assumption, the injected 
electron flux is downward beamed (/io = 1), and its en- 
ergy spectrum (see Fig. ^ is a broken power law oc E~^ 
with an index (5 = 6 above and 5 = —2 below the cutoff 
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TABLE 1 

Summary of simulation runs. 

tluns Injected Electron Particle i-max ^umax ^min tD>ioo tapcx Tmax rie.max 

Spectrum Ang. Distr. Transport (kms~^) (s) (kms~^) (s) (s) (10^ K) (10^° cm~^) 
0(ld) power law beamed analy. approx. 565 35 —115 10 29 21.1 6.96 
H(ybrid) power law beamed Fokker-Planck 570 35 -106 9 28 22.0 7.44 
N(ew) stoch. accel. isotropic Fokker-Planck 598 36 -102 9 23 25.9 8.44 

Notes. Runs O and H have an identical injected power-law electron flux, with a spectral index <5 = 6 and low-energy cutoff Ec = 15 keV; 
Kmax and ty^^^: maximum (upflow, -u > 0) velocity and its time stamp; Vmin: minimum (downflow, < 0) velocity (in the upper 
chromosphere); ti,>ioo- time when the upflow velocity exceeds lOOkms"^; fapcx: time when the evaporation front (density jump) reaches 
the loop apex; Tmax and rie^max: maximum coronal temperature and electron density. All runs have the same peak energy deposition 
(electron heating) flux for the loop as a whole, -Fmax = 2.67 X 10^" ergss"^ cm~^. 




1 10 100 1000 10000 

E (keV) 



Fig. 4. — Electron flux spectra F(E) times at the peak time 
{t = 30 s): (1) acceleration region flux Fac, escaping flux Fosc, and 
equivalent thick-target flux -Fthick for Run N, and (2) injected flux 
for Runs O and H ((5 = 6 above Ec = 15 keV). The triple-dot- 
dashed line (arbitrary scale) represents the source term Q{E) of a 
Maxwellian (thermal, fcjT = 1.53 keV) distribution used in Run N, 
which peaks a.t E = 3fci,T in this E^F(E) plot. 

("knee") energy = 15 keV. Here the differences are: 
(1) our peak energy flux, i.e., parameter F in equation (9) 
of MEL89, is 2.67 x 10^° ergsem-^ s-\ while they used 
5 X lO^*' ergsem"^ s~^, and (2) we assume a fully ionized 
hydrogen plasma while they included helium of 6.3% of 
hydrogen in number density. The latter difference only 
changes the absolute mass density by 11%, while the rela- 
tive differences between our models may not be affected. 

4.2. Run H: Hybrid Model 

Here we used our Fokker-Planck transport code in 
place of the approximate analytical expression used 
above to evaluate the heating rate. We injected elec- 
trons of an identical power-law spectrum with a narrow- 
Gaussian (cTp = 0.01) pitch-angle distribution to emu- 
late the beamed distribution in Run O. The main differ- 
ence here is that the transport code properly treats the 
diffusion process of pitch angle change due to Coulomb 
collision. For the particle code, the energy space is di- 
vided into 200 uniform logarithmic bins in the range of 
511 X [10"'^, 10^] keV, while t he pitch a ngle space is di- 
vided into 100 (vs. 24 used in lLiull2006[ ) uniform bins in 
the [0, tt] range. 

4.3. Run N: New Model 

This is a typical simulation using our new model. It 
is the same as Run H, except that the injected beamed 
power-law electron flux is replaced with an isotropic flux 



given by the stochastic acceleration code. We used the 
same acceleration parameters as PL04 (see their Fig. 12), 
i.e., the characteristic acceleration rate = 70 s~^, 
Ue = 1.5 X 10i°cni-3^ Bo = 400 G, hT = 1.53 keV, and 
the acceleration region size L = 5 x 10® cm. We modu- 
lated the rate {Q, see eq. [T]) of electrons being supplied 
to the acceleration region with the same triangular time 
profile, such that the peak electron heating flux i^max 
equals that of Run O or H. 

Figure |4] shows various electron flux spectra used in 
this study. In comparison with the background thermal 
distribution, both the acceleration region flux i^ac and 
escaping flux i^osc have a quasi-thermal component that 
smoothly extends to a nonthermal tail at high energies, 
i^osc is smaller than i^ac and their relative difference de- 
creases with energy due to the energy-dependent confine- 
ment of electrons by turbulence in the acceleration region 
(see eq. [5]). Unlike that (dot-dashed) in Run O or H, the 
flux i^osc injected into the transport region does not in- 
voke any arbitrary low-energy cutoff. The two fluxes, 
however, have similar slopes in the intermediate energy 
range around 20 keV. The equivalent thick-target elec- 
tron flux i^thick (see as expected, has a harder spec- 
trum than Fac and i^osc in the 10-1000 keV range. 

5. SIMULATION RESULTS: COMPARISON OF FLUID 
DYNAMICS 

To determine how much our proper transport and ac- 
celeration calculations affect the atmospheric response, 
we compare the results of the three simulations, using 
Run H as the reference case. 

5.1. Hydrodynamic Evolution 

As shown in Figure [51 Run H exhibits similar general 
HD evolution as described in MEL89 (their Figs. 1 and 
2). Electron heating (5e) is initially concentrated in the 
upper chromosphere, producing overpressure that drives 
an upflow {v > 0) and a recoil downflow {v < 0). At 
t = 9s, the upflow velocity exceeds lOOkms"^ and a 
density jump or evaporation front has developed slightly 
above the TR. It travels upward and reaches the loop 
apex at i = 28 s. The density jump is then reflected back 
and the material piles up due to the reflective boundary 
condition imposed, which can be understood as plasma 
flow from the other half of the loop in the assumed sym- 
metric geometry. The upflow reaches its maximum veloc- 
ity of Umax = 570kms~^ at t = 35 s, which is delayed by 
5 s from the energy deposition peak at < = 30 s. Chro- 
mospheric evaporation then gradually subsides. These 
features of the temporal evolution can also be seen from 
the history of various quantities at a fixed position in 
the upper corona as shown in Figure [6l Note that, late 
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Fig. 5. — Evolution of electron density, temperature, pressure, 
velocity, and electron heating rate for Runs H and N. 

during the simulation, the coronal temperature gradu- 
ally decreases mainly due to conductive cooling, while 
the coronal density continues to increase, even after the 
cease of electron heating at i = 60s. This is caused by 
sustained chromosphcric evaporation that results from 
heating of the chromosphere by the same conductive flux 
that cools the hot corona. 




40 60 
Time (s) 

Fig. 6. — History of electron density, temperature and velocity 
at s = 1 Mm from the injection site for different simulations. 



In comparison, we find that Run O is very similar to 
Run H (see Fig. [6] and Table [T]) but less intense. This 
can be more clearly seen from various HD quantities and 
heating and cooling rates during the early phase shown 
in Figure [T] The small differences (on the order of 10%) 
between the two cases are mainly caused by the slight 
overestimate of electron heating of Run O in the chromo- 
sphere (Fig.[71i), due to its inaccurate way of calculating 
particle transport noted earlier. This indicates that, for 



this specific case, the analytical heating function used 
in MEL89 provides an acceptable approximation for the 
more accurate Fokker-Planck calculation. 

In contrast, the HD evolution of Run N is faster and 
more intense than that of Run H (Figs. [5] and [H Table[T]), 
despite the same maximum electron heating rate for the 
two cases. These differences are primarily due to the 
different energy spectra of injected electrons, further 
enhanced by their different pitch-angle distributions: 

(1) In Run N, electrons of a few kcV in the quasi- 
thermal component of the spectrum (Fig. [4]) are small in 
energy but large in number and thus dominate the total 
energy content. These electrons produce heating at rela- 
tively small column depths by Coulomb collisions. This 
occurs high in the corona (Fig. [Tiff); where the radiative 
loss rate irad (Fig- IZls) is relatively small due to the low 
density and high temperature. As a result, significant 
net heating sets in there, which leads to a local temper- 
ature (Fig. 13)) and pressure surge. This local coronal 
heating is enhanced by the large effective column depth 
(iVoff = N/ (fi) , where (/i) is the mean pitch angle cosine) 
resulting from the isotropic pitch-angle distribution of 
the injected electrons. The increased temperature leads 
to a large downward heat conduction flux (Fig. [Tlf), and 
the pressure gradient force drives a downward mass flow 
in the high corona (Fig.[7]c). 

(2) In Run H, contrastingly, the electron spectrum 
(Fig. S]) peaks at the cutoff energy Ec, which leads to 
a deficit in low-energy electrons. In addition, the pitch- 
angle distribution here is beamed (rather than isotropic). 
This electron population, on average, penetrates deeper 
into the atmosphere than that in Run N and deposit 
their energy primarily in the upper chromosphere. This 
results in less heating in the corona and stronger and 
more widespread heating in the chromosphere (Figs. [3? 
and[71i). Consequently, in spite of the larger and broader 
radiative cooling (Fig. [7]e) , the local overpressure in the 
chromosphere is stronger than that in Run N early on, 
which drives a higher velocity upflow (Fig. [Tjc, t = 6 s). 
Also, unlike Run N, there is no significant downward 
coronal heat conduction (Fig. [3f) or mass flow. 

5.2. Heating and Cooling 

A remaining question is why Run N has more dramatic 
overall HD changes in the long run. To answer this, we 
examine the relationship between different energy gain 
and loss terms — electron heating, radiative loss, and 
conductive heating and cooling — particularly early in 
the flare and near the TR where chromosphcric evapora- 
tion takes place. 

(1) In Run N at i = Is (see top panel of Fig. [8h), 
the electron heating rate Se peaks in the TR because 
of the sharp increase there in ambient density and asso- 
ciated collisional energy loss of energetic electrons. So 
does the radiative loss rate Lrad, since it is proportional 
to UeUp and (f>(T) that peaks at T ~ 1 - 3 x 10^ K which 
is the TR temperature. However, due to their differ- 
ent functional dependencies on density and temperature, 
Se peaks at a slightly lower position than Li-ad- Their 
combination Se — ./^rad {panel 2, Fig. ^p) thus results in 
cooling in the upper TR and heating in a shallow layer 
below it in the upper chromosphere. Meanwhile, the con- 
ductive flux carries energy from the hot upper corona 
to the upper TR where localized heating (S'cond) is pro- 



Combined Modeling of Acceleration. Transport, and Hydrodynamics in Solar Flares 



9 




r (d) Electron Heating 
















- I 


\ : 







4 6 
s (Mm) 




4 6 
s (Mm) 



Fig. 7. — Same as Fig. [6] but for snapshots of various quantities as a function of distance at selected times early in the flare. 



duced and counteracts radiative cooling. (Conduction is 
prohibited in the chromosphere where the temperature is 
maintained at 10^ K.) The net energy gain resulting from 
the interplay of electron heating, radiative cooling, and 
heat conduction is thus localized in the upper chromo- 
sphere, where temperature is raised substantially (panel 
3, Fig. (SJa). This leads to a local pressure hump which 
drives an upflow into the corona and a downflow into the 
chromosphere (panel 4, Fig.lSh). 

As time proceeds (see Fig. \8i}) and chromospheric ma- 
terial is being heated from T = 10* K to ~10^K where 
$(T) reaches its maximum, radiative loss gradually over- 
takes electron heating in the TR and upper chromo- 
sphere. This means that energy directly deposited by 
electrons in these places is immediately radiated away 
and very little is left to heat the plasma. However, as we 
noted above, a significant portion of the energy content 
of the injected electrons in Run N is deposited in the up- 
per corona (Fig. [TH) where radiative loss is negligible and 
then transported by heat conduction to the lower atmo- 
sphere. In time, conduction plays an increasingly impor- 
tant role in heating the lower corona and TR as it gradu- 
ally exceeds the direct net heating or combined electron 
heating and radiative loss Se — irad (Fig- E?') in these 
regions. Because of this, as of i = 7s the location of the 
primary net heating, i.e., the highest local overpressure, 
and the maximum downflow velocity have shifted from 
the upper chromosphere up to the TR. In fact, the peak 
of the conductive flux i^cond (see Fig. ^f) and the region 
of significant conductive heating (S'cond oc —dFcond/ds) 



below it are initially located in the upper corona and 
propagate down the loop with time. As the i^cond peak 
approaches the TR at ^10 Mm during the interval of 
8-10 s, the maximum upflow velocity Wmax in the loop 
increases abruptly (see Fig. [9|) . This further emphasizes 
the role played by heat conduction here in redistribut- 
ing energy deposited by electrons and in driving chromo- 
spheric evaporation. 

(2) In Run H, the injected electron flux with a low- 
energy cutoff has profound consequences. As noted 
earlier, lack of low-energy electrons makes the chromo- 
sphere, rather than the upper corona, the primary loca- 
tion of direct heating early in the flare (see Fig. [TH). 
Net direct heating 5*6 — Lrad and the increase of lo- 
cal temperature and pressure extend from the TR to 
much deeper layers of the chromosphere than in Run N 
(Fig. [He). Moreover, since the coronal temperature does 
not increase rapidly, the downward conductive flux here 
is more than an order of magnitude smaller (Fig. [30 ■ 
Net direct heating generally dominates over conductive 
heating when integrating over the volume of the lower 
atmosphere. As a result, in comparison with Run N, a 
relatively larger portion of the total energy content of the 
injected electrons is lost in radiative cooling. This is why 
the overall HD development here is more gentle than that 
of Run N, despite the fact that they have the same en- 
ergy deposition flux. The primary underlying physics is 
their different spatial distributions of electron heating Se 
caused by their different electron injections. Note that 
MEL89 obtained qualitatively similar results when dif- 
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(a) Run N (t=1 s) 



(c) Run H (t=1 s) 
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Fig. 8. — Comparison between Runs N (left) and O (right) of a 
detailed view near the TR at two selected times. At each time, the 
quantities in the four panels (from top to bottom) are: (1) electron 
heating rate Se and radiative loss rate Lrad, (2) Se — irad (dotted) 
and conductive heating rate Scond (solid), (3) electron density (left 
scale) and temperature (right scale), and (4) pressure (left) and 
velocity (right). Note different vertical scales. 




Fig. 9. — Early history of the maximum upflow velocity I'max 
(left scale) and the position of the maximum conductive flux 
(right scale) for Run N. The two vertical lines mark the interval 
when I'max experiences a sharp increase while the maximum Fp^nd 
rapidly approaches the TR. 



fcrent values of the cutoff energy or spectral index were 
considered. We note that, later during the impulsive 
phase in Run H (Fig. (5]), as the coronal density has in- 
creased considerably due to chromosphcric evaporation, 
relatively more electron energy is directly dumped in the 
corona while less in the chromosphere. Consequently, 
conductive heating in the TR becomes important. In 
this sense, the physical distinction between the two runs 
gradually diminishes in the late stage of the flare. 

5.3. Velocity Distribution 

Here we examine observables that can be checked 
against data: (1) the temperature dependence of the 
plasma velocity (Fig. [TOl left), and (2) the velocity dif- 
ferential emission measure (V DEM; Fig. [TOl right) de- 
fined by iNewton et all (|1995( ) as the emissivity {G[T]) 
weighted measure of material with line-of-sight (LOS) ve- 
locity vi^os ■ Assuming that the flare is located at the so- 
lar disk center such that the LOS is perpendicular to the 
loop apex, we calculated the specific VDEM/a(s) for the 
Ca XIX 3.18 A hue following [Newton et al.lHere GfT) is 



given by the Chianti package (jYoung et al 
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Fig. 10. — Comparison between Runs H and N: velocity vs. tem- 
perature { left) and corresponding specific velocity differential emis- 
sion measure (photons cm~^ sr~^ [km s~^]~^) vs. LOS veloc- 
ity (right). 



Early in the flare (Figs. [TOb and [TOb). Run H has 
higher upflow velocities and downflow temperatures than 
Run N, owing to its stronger net direct heating with a 
deeper extent in the chromosphere (Fig. \8jp). Run N 
has an additional hot downflow component at >2 MK. 
due to the expansion of the heated upper corona men- 
tioned earlier. This hot component also dominates the 
VDEM at wlos < (Fig. [TOb ) because of the sharp rise 
of G(T) up to T = 29 MK. Later at t = 28 s (Fig. [TOF*), 
Run N overtakes H in upflow velocity and temperate, 
while its hot coronal downflow has disappeared. Run N 
has a bimodal VDEM (Fig. fTOH ) with a strong station- 
ary, hot component located near the loop apex (also see 
Fig. [5]), while Run H exhibits a more gradual progression 
toward high velocitie s. This distinction is similar to what 
iNewton et al.l (|1995l . see their Fig. 2) found for models 
with different initial coronal densities. 
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Fig. 11. — Spectra of angle-integrated electron flux (top) and corresponding photon emission rate {bottom) multiplied by energy squared 
at three selected times for Run N. The thick dashed line represents the LT (acceleration region) spectrum (Fac or /lt), while the thin 
colored lines [solid, dotted, dashed, dot-dashed, and triple- dot-dashed) are the spectra at distances s = 4, 10, 11, 12, and 13 Mm from the 
injection site. The thick dotted line indicates the escaping electron flux (Fcsc) in the top panels and the equivalent FP photon emission 
rate (Ipp) in the bottom panels. The legends show the current values of the position (str) and column depth (Ntr) of the TR. 



6. 



EFFECTS OF HYDRODYNAMICS ON PARTICLE 
TRANSPORT AND X-RAY EMISSION 



We now turn our attention to the effects of fluid 
dynamics on particle characteristics, namely, electron 
transport and nonthermal radiation. Here we present 
the result of Run N as a typical example. We will 
examine first the energy distributions of electrons and 
bremsstrahlung photons, and then their spatial distribu- 
tions. 

6.1. Electron and Photon Energy Spectra 

Figure [Tl] (top panels) shows the evolution of the angle- 
integrated electron flux spectrum E'^F{E, s) at different 
locations in the loop. In general, at a given time and 
a moderately large distance or column depth, there is 
a deficit in low-energy electrons due to coUisional en- 
ergy losses and scatterings on the paths of the electrons 
from the injection site. This appears as a turnover in 
the spectrum and slight spectral hardening just above 
the turnover energy. As distance increases, progressively 
more low-energy electrons are lost, and thus the overall 
flux decr eases while the spectral tu rnover shifts to higher 
energies (|Leach fc Petrosian|[T98l[) . At t = 2s (Fig. [Tib) 
the TR is located at distance Str = 9.97 Mm and column 
depth A^^tr = 6.1 x 10^^ cm"^. The two fluxes at s = 4 
and 10 Mm [thin solid, dotted) are located in the low- 
density corona or TR at small column depths from the 
injection site, and thus appear similar in shape to the 
escaping flux Fcsc {thick dotted). Other fluxes at s = 11, 
12, and 13 Mm are located in the chromosphere at large 
column depths, and thus exhibit substantial reduction of 
low-energy electrons. 

In this study, the flux (-Fcsc) injected into the transport 



region does not change with time in spectral shape but 
only varies in normalization. So does the electron flux 
at a given column depth. However, as chromospheric 
evaporation develops, the TR retreats to lower altitudes 
(Fig. [5]), while the coronal density increases. This causes 
the change of the column depth and thus the electron 
spectrum at a fixed location, which is what we notice 
here. At i = 30 s (Fig. when the TR shifts down to 
Str = 11.3 Mm at Ntr = 4.2 x lO^^cm-^ the spectrum at 
s = 1 1 Mm looks similar to the other coronal spectra (at 
s = 4 and 10 Mm) but different from the chromospheric 
spectra. Meanwhile, the relative difference between the 
first coronal spectrum (s = 4 Mm) and the injected spec- 
trum i^csc becomes larger, because of the increasing coro- 
nal density and column depth between them. This trend 
continues through the end of the simulation. 

The bottom panels of Figure [Tl] show the correspond- 
ing bremsstrahlung photon spectra e^/(e, s) defined in 
equation ([7]). Like its electron counterpart, the LT pho- 
ton spectrum e^/lt {thick dashed) shows a thermal-like 
component at low energies and a nonthermal tail at in- 
termediate energies with a cutoff at high energies. The 
equivalent FP spectrum e^/pp below str {thick dotted, 
defined in i i2.4|) is harder than the LT spectrum. Their 
shapes in the lOs-lOOs keV range are commonly ob- 
served for LT and FP sources. At larger distances, 
the spectra {thin lines) become progressively harder be- 
cause the corresponding electron sp ectra have the same 
variation (jLeach fc Petrosiaiilll983[ ) as we have noticed 
above, and high-energy (>10s keV) emission increases 
because of higher ambient densities. As the coronal den- 
sity increases with time, for the same reason mentioned 
above photon spectra at positions originally located in 
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Nj^=4.2e+19 cm'' 
' Sj^=1 1.3 Mm 




e = 3.1 6.1 12.3 24.5 
48.8 97.4 294.1 keV t 



5 

s (Mm) 

Fig. 12. — Spatial distributions of angle-integrated electron flux {top) and photon emission rate [bottom) at three selected times for 
Run N. From top to bottom, different lines in each panel represent electron or photon energies of 3.1, 6.1, 12.3, 24.5, 48.8, 97.4, and 
294.1 keV. The step in the s <0 region corresponds to one half length of the acceleration region at the top of the loop (see Fig. [TJ. The 
step at the TR is due to the jump in the ambient density. 



the corona become progressively harder. At the same 
time, as the TR treats, spectra at positions originally in 
the chromosphere but now in the corona become softer. 
In other words, the difference in spectral shape between 
X-ray emissions at different positions diminishes with 
time (Figs.HHs andiny). 

6.2. Electron and Photon Spatial Distributions 

Figurc fT^ f^op panels) shows the evolution of the angle- 
integrated electron flux F{E, s) as a function of distance 
at different energies from 3 to 300 keV. The step at s = 
is owing to the difference between the LT (acceleration 
region) flux i^ac and the escaping flux i^esc mentioned 
before (see Fig. |4]). In general, the electron flux de- 
creases with distance or column depth from the injec- 
tion site. The slope of the curve -~dF{E, s)/ds is steeper 
at lower energies because lower energy electrons lose en- 
ergy faster and are more sensitive to pitch-angle scatter- 
ing. At a given energy, the slope depends on the ambi- 
ent density, because dF{E, s)/ds = ne[dF{E, N)/dN], 
where F(E, N) is generally a smooth function of N 
(jMcTiernan fc PetrosianI 1 1990( 1. Therefore, if there is 
a rapid increase in density with distance, the slope in- 
creases quickly; the opposite would be true if the density 
were to decrease sharply. A sharp break is obvious at the 
TR, and a milder break occurs at the evaporation front 
(e.g., at s ~ 5 Mm in Fig.[T2b). An upward break or flat- 
tening is visible at s ~ 3 Mm where the density jump is 
reflected back from the loop apex at t = 30 s (Fig. [T2fc ) . 
It is also visible just below the TR (see Figs. [TTb and 
112b ). where a sharp density decrease from the narrow 
density spike occurs (Fig. [5]). As density increases due 
to chromospheric evaporation, the spatial distribution in 
the corona becomes steeper. The slope variation with 



distance (i.e., breaks) and time is more pronounced at 
lower energies for the same reason noted above. 

Figure [12] ( bottom panels) shows the evolution of the 
corresponding spatial distribution of angle-integrated 
bremsstrahlung emission /(e, s). In the early stage, low- 
energy emission comes primarily from the LT, while high- 
energy emission is concentrated below the TR. Because 
nonthermal bremsstrahlung radiation is proportional to 
both the electron flux and local target density, photon 
emission reflects details of the density distribution in a 
more pronounced way than the electron flux profile. As 
is evident, the emission profile clearly indicates density 
features, including the evaporation front and the density 
spike just below the TR. The local emission enhancement 
at the evaporation front moves upward with time, which 
may be responsibl e for the observed X-ray sources mov- 
ing along the loop (|Liu et al.ll2006bl : ISui et "aI1l2006[ ) men- 
tioned in § 11.11 As time proceeds, relatively more emis- 
sion comes from the coronal portion of the loop because 
of the increased density. Specifically, at low energies, the 
emission intensity decreases with distance in the corona 
more sharply than before. At intermediate energies, we 
find a temporal transition from FP-dominated emission 
to LT-dominated emission, which occurs at progressively 
hi gher energies, remi niscent of the phenomenon observed 
bv lLiu et all (|2006b[ ). At high energies, such a change is 
not visible because the high density corona still looks 
transparent to high-energy electrons, but the retreat of 
the TR is obvious. 

7. SUMMARY AND DISCUSSION 

We have performed simulations of solar fiares that self- 
consistently combine acceleration, transport, and radi- 
ation of energetic electrons (using the Stanford unified 
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code) with fluid dynamics of the atmospheric response 
(using the NRL flux tube code). As the first success- 
ful one of its kind, this model improves on previous HD 
simulations in two major aspects. First, it includes more 
accurate evaluation of electron heating from full Fokker- 
Planck calculation of particle transport. Second, it uses 
a more realistic electron spectrum from the stochastic 
acceleration model (PL04) as the injection to the trans- 
port region. We compare this more advanced treatment 
with models in which an ad hoc electron distribution of 
a power-law with a low-energy cutoff is injected into the 
loop and/or transport is dealt with approximately. Our 
conclusions are: 

1. For the specific injection of beamed, power-law elec- 
trons, the old analytical model of MEL89 (Run O) pro- 
vides an acceptable approximation. Its result differs by 
~10% from that of the reference hybrid model (Run H) 
obtained by the more accurate Fokker-Planck calculation 
(see Table [D and Figs. [5] and [T]) . 

2. In the new model (Run N), where the injected elec- 
tron spectrum is based on stochastic acceleration, we find 
higher coronal temperatures and densities, larger upflow 
velocities, and faster increases of these quantities than 
the hybrid model (Run H, Fig.[S]). This is mainly because 
the new injected electron spectrum smoothly spans from 
a quasi-thermal component to a nonthermal tail (Fig.|3]). 
The low-energy electrons in the quasi-thermal regime, 
which contain the bulk of the total energy budget, de- 
posit their energy primarily in the corona. This results 
in significant coronal heating and thus a large downward 
heat conduction flux that helps drive "evaporation" of 
plasma at the TR. In contrast, the electron spectrum in 
the hybrid model with a low-energy cutoff leads to more 
energy directly deposited in the chromosphere, which is 
radiated away more quickly, leaving less energy to pro- 
duce actual heating (Figs. [7] and [8]). This is qualitatively 
consistent with the conclusion of MEL89 where an elec- 
tron spectrum with a smaller low-energy cutoff or steeper 
slope resulted in a stronger chromospheric evaporation. 

3. The energy and spatial distributions of energetic 
electrons and bremsstrahlung photons bear the finger- 
print of the changing density distribution caused by chro- 
mospheric evaporation. In general, as time proceeds, 
the electron and photon spectra at positions remaining 
in the corona become progressively harder because of 
the increasing coronal density, while those at positions 
previously in the chromosphere and now in the corona 
(due to the retreat of the TR) become softer (Fig. [TT|) . 
Any density jump in space results in a sudden change in 
the spatial distributions of energetic electrons and X-ray 
photons (Fig. I12p . In particular, the evaporation front 
appears as a local emission enhancement, which, in prin- 
ciple, can be imaged by X-ray telescopes. 

7.1. Comparison with Observations 

Ov er several decades (see review by lAntonucci et al.l 
|1999() . Doppler observa tions have indicate d hot, fast 
(>100kms"^) upflows (jPoschek et al] 119801 an d cool- 
slow (>10kms~^ ) dow nflows (jWuelser et all 119941 : 
iBrosius fc Phillips! |2004[ ) during flares. This IS con- 
sistent with chromospheric evaporation and momen- 
tum recoil shown i n our and earlier si mulat ions (e.g., 
Mariska et all 119821: iFisher et al.l Il985c[ ). In addition^ 




multiple temperatures obtained from Hinode EUV Imag- 
ing Spectrometer (EIS) observations, which show ex- 
cellent agreement with the Run N curve in Figure llOb 
throughout the 10^-10^ K range. Such an agreement 
with HD simulations has not been seen before. 

Heat conduction, compared with more popular direct 
electron heating, plays an important role in driving chro- 
mospheric evaporation in our new mode l. Observational 
suppo rt of this was first reported by iZarro fc LemenI 
()1988f ). who found an upflow velocity < 50k ms~^ at 
_6MK during the cooling phase of a flare. iMilliganI 
recently reported an unusually high temperature 
(2 MK) downflow at ~14kms~^ in a microflare with 
no detection of HXRs (implying a low flux of nonther- 
mal electrons), which possibly results from conductive 
heating. This downflow could be due to the thermal 
expansion early in the corona (Fig. 1 10b ) or the momen- 
tum rec oil later in the chromo sphere (Fig. \Wh). More 
recently, iBattaglia et al.l ()2009[ ) interpreted the growths 
of the SXR emission measure observed early during 
flares (before HXRs being observed) as new evidence of 
conduction-driven chromospheric evaporation. 

Our simulations predict that evaporation upflows tend 
to have higher temperatures when conductive heating 
dominates over direct electron heating, while the oppo- 
site is true for recoil downflows (Fig. fTU]) . This can be 
checked against observed distri butions of the plasma bul k 
velocity vs. temperature fe.g.. iMiUigan fc Dennis! [2009[ ). 
Our simulated VD EM can also be readil y used to synthe- 
size emission lines (jNewton et al.|[l995D to be compared 
with observations and help differentiate theoretical mod- 
els, because of the sensitive dependence of VDEM on 
heating mechanisms. 

7.2. Future Work 

This paper is the first in a series and we have es- 
tablished the numerical model and presented initial re- 
sults. In foUowup papers, we will explore the parameter 
space and use this model to investigate t he Neupert effect 
and the observed moving X-ray sources (iLiu et al.ll20"06bl : 
ISui et a l. 2006). More importantly, we will incorporate 
the atmospheric feedback on the acceleration process. 
This is because chromospheric evaporation may change 
the physical condition (e.g., plasma density and tempera- 
ture) in the LT acceleration region. The density enhance- 
ment, for example, causes the ratio of electron plasma 

1 /2 

frequency to gyro-frequency a = ujpc/fle Ue / Bq to 
increase. This can lead to the reduction of the efficiency 
of electron acceleration (PL 04) and thus the quench ing 
or spectral softening (e.g., iParks fc Winckleil Il969f ) of 
nonthermal HXR tails observed during the late stages 
of flares. 

Some technical aspects of this model can be improved 
in the future: (1) The fully ionized hydrogen plasma 
assumed here can be replaced by a plasma of a more 
realistic solar abundance, with the inclusion of neu- 
trals and ionization equilibrium. (2) The "cold" tar- 
get assumption in the transport code can be abandoned 
and replaced with Couloinb diff usion in energy space 
(|SpitzeH[l96l iMillcr et al." 19960 for a general "warm" 
target plasma (jEmslic. 2003.) . (3) In particle transport^ 
and HD calculations, one can include energetic protons 
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(c.f., lEmslie et al.l |1998[ ) and heavier ions, whose mo- 
mentum loss to the background p lasma, in addition to 
the overpressure (|Kosovichevll2006[ ) produced by electron 
heating, could contribut e to generating seismic waves 
observed in some flare s (jKosovichev fc Zharkoval Il998t 
iDonea fc LindsevI 120051 ). (4) In the long run, we intend 
to implement time-dependent particle transport calcu- 
lation, full loop simulation in an asymmetric geometry, 
return currents, and radiative transfer. 

The combined treatment of the particle and fluid de- 
scriptions of plasma presented here opens a door to a 
broad range of applications. This model, with proper 
modifications, can be applied to environments where in- 
terrelated particle acceleration and transport and plasma 
flows are present, such as (exo)planetary auroras (e.g.. 
iLiu fc Airapetianll2008f ) and flares on other stars and in 
accretion disks near black holes. 
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